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Two Sample Covariances from a Trivariate Normal Distribution 
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Abstract. The joint distribution of two ofF-diagonal Wishart matrix 
elements was useful in recent work on geometric probability [12j. Not finding 
such formulas in the literature, we report these here. 

Let {A, B, C) be trivariate normally distributed with known mean (0, 0, 0) and 
covariance matrix 



Gov 





1 - ^2 > 0, 



1 - 2p2 + a > 0. 



We wish to evaluate the probability that AC < x and BC < y. On the one hand, 
the joint density of {AC, BC) is pL] 



f{x,y) 



1 '"" Vf 



- (px + py- v^\/(l - p2)a;2 - 2(-p2 + a)xy + (1 - p'^)y^) 



2n 



1 — p2)x2 — 2(— p2 + a)xy + (1 — p^)y'^ 



where ^ = l-2p2+a, r] = (l+a)/(l-a). On the other hand, {A^, 2AB, 2AC, B\ 2BC, C^) 
is pseudo- Wishart distributed with 1 degree of freedom and thus {AC, BC) has char- 
acteristic function [2] 



F{u,v) 



det 




/ 1 



2i 



a p 



alp 

1 



u/2 
v/2 
u/2 v/2 



-1/2 



1 + (1 — p2)n2 — 2ipu + (1 — p2)^2 _ 2ipv + 2(— p^ + (j)uv 

Our first task is to confirm that f{x,y) and F{u,v) are indeed a Fourier transform 
pair. 
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Now let {Ai,Bi,Ci), (A2, -B2, C'2), • • •, {An,Bn,Cn) be a random sample from 
A^(0, S) and define sample covariances 



7A,c = Yl ^iQ' 7b,c = Y. Bj^: 



The marginal density for ^a,c is well-known [H IH O El O [H [H [TO] : 

^.^l-p'V(nl2r ^ ir^^j ^'"-')/^ (737 j 

and likewise for 7b,c, where i^^^_jj^2(^) is the modified Bessel function of the second 
kind. The joint density fn{x,y), however, is not fully understood [11] even though 

Fn(u,v) = j7^ 

[1 + (1 - p^)u^ - 2ipu + (1 - p2)t,2 _ 2ipv + 2(-p2 + a)uvr^ 

is comparatively simple. We shall determine fn{x,y) for n = 2,3,4 and then for 
arbitrary n, closing an evidently open issue. A special case involving f^{x,y) was 
examined in |12j, for which p = a = 1/2, to answer a geometric probability question. 
Our discussion generalizes this earlier work. 

1. Fourier Transform Pair 
Our objective is to evaluate the integral: 

00 00 exp 7 ((p + i^u)x + (p + i^v)y — ^y/ax'^ — 2bxy + ay'A 



^{u,v) = — / / ^ / o , , ^ '-dxdy 

In J J y/ax'^ — bxy + ay^ 



-00 —00 



where a = 1 — p^, 6 = — p^ + a. Let 

r T 

X = —r= {\cos{6) — Ksin(6')) , y = —= (Acos(6') + fi;sin(^)) 
v2 v2 

then ax^ — bxy + by'^ = r^ under the requirement that 

A ' 



The Jacobian determinant is 



1 V 



,V2. 



A cos{6) — K sin{6) r (—A sin(0) — k, cos{6)) 
A cos(6') + K sin{9) r (-A sin(6') + k cos{9)) 



Xnr 
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hence dx dy = Xnr dr d9. We obtain 

27r oo , 

\? V 



+ (p + i^v) (A cos(6') + K. sin(6')) - J2r] 



dr dO 



27r oo / 

Xk f f I 1 r 







A< 



2n 



2\p cos{9) + iA.^('U + v) cos{9) — iK^iu — v) sin(^) 
1 



dr dO 



y/2Ti J 2Xpcos{6) + iX^{u + v) cos{6) - iK^{u - v) sm{6) - ^fh] 



dd 



since 



Re 



2Xpcos{6) + iX^{u + v) cos{6) — iK^{u — v) sin(6') — J2ri 



= 2Xpcos{e)- ^2ri 
< 2Xp - ^/2r] 



and 

,oo 1 o 1 + cr 2 

4AV -2v = 4- / - 2- = 

1 — 0" 1 — a I — a 

Let z = exp(i 9), then d9 = —i dz/z and 
iXnC, 



l-2p' + a) < 0. 



F u,v) 



dz 



v/27r ^c 2Api(2; + i) + ^A^l^ + t;)|(2 + i) - t^^iu - v)^{z - )) - ^m z 



^iXni 



TT Jc 2Ap(z2 + 1) + iX^{u + t;)(22 + i) _ /.^(^ _ ^)(^2 _ i) _ 2^z 
The denominator of the integrand can be rewritten as 

(2Ap — K^u + K^v + iX^u + iX^v)z^ — 2J2riz + (2Ap + k^m — k^v + iA.^-u + iX^v). 
Let (5 = 2p/^. The two poles Zpos, -Zncg of the integrand are 



dz. 



2y/2ri ± ^J8r] - 4(2Ap - k^u + <t; + iX^u + iA^w)(2Ap + <m - k^v + iX^u + iX^v) 

2(2Ap — K^M + K^v + iX^u + iX^v) 



2r] ± ^A2^2(m + t; - i5)2 + k'^^^{-u + v)^ + 2r] 
2Xp — K^u + K^v + iA^-u + iA.^f 
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and Zj^eg is always inside C, Zpos is always outside. Clearly z^^g is a pole of order 1 
and the associated residue is 

lim 



Z->2neg (2Ap — K^U + K^V + iX^U + iX^v){z — Zpos) 

1 1 



2 ^A2^2(m + V- i5y + /€2^2(-M + v)^ + 2r] 

multiplying by (2ni){\/2i\K^/7i) gives 

V2\k^ 

V2Xk 



JX^iu + V — i5y + k'^{—u + vY + a^ 

where a = ^/2ri/C,■ This is the most useful expression for our purposes. The original 
formula for F{u, v) can now be confirmed. 

For example, if p = a = 1/2, we have ^ = 1, ?7 = 3, A = -\/2, k = 1, 6 = 1, 
a = \/6. In this special case, our expression for F(u, v) becomes 

2 1 



2{u + v- iy + {-u + vy + 6 J(f - i)(3f + i) + (I - 2)(3f + i) + 2f | - 1 



and this is consistent with [T2] . 

We can also argue in reverse, that is, evaluating instead the integral: 

J — OO CXD 

V2Xk f f exp{-iux - ivy) 

j{x,y) = - — r— / / ^av du. 

(^"^^ -OO -OO \/^^(« + ^ - ^^Y + >^^{-u + ^Y + "^ 

Let s = u + V, t = —u + f , then u = {s — i)/2, f = (s + t)/2 and 

s t s t 

—lux — ivy = —i-x + i-x — i-y — i-y 
^ 2 2 2 2 

X + y —X + y 

= —I s — I 1: 

2 2 

hence 

fU y) = ^^ I I ^S^l=^L^^jLdtds 

'' ''' V2{2iTy_l_l ^x'^s-iSy + KH' + a^ 
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since the Jacobian determinant is 1/2. The inner integral is [13 

—X + y 



exp(-i^t 



2^^ /I 
^dt = — Ao — 



\^{s — iby + d' 



Define 



then the outer integral is [13 



1 

7 = - 



—X + y 



oo 

/ Kq {-i^\^{s-ibY^o?\ exp i-i 



X + y 



s I c/s 



'x + y 
IT exp I — r — 



exp -^.yA2 + ^ 



ivTKexp 



.p ^-^J^^x + yf + XHx-yf^ 
2 "^ Y^/s:2 {x + yf + \^ {x - yf 

Multiplying by {\k / {^/2{2■Kf)) ■ (2/k), we obtain 



V^' 



exp 



vr 



x + y 



^^v[-^\Ji^'^{x + y) +X^{x-yy 

\J K? {x + y) + A^ {x — y) 



and the original formula for f{x,y) can now be confirmed. For example, if p = a 
1/2, the probability that both AC > and 5C > is 1/2 = 0.5. 



2. Sample Size n = 2 



We wish to evaluate 



f2{x,y) 



2,^2 



exp [-i^s - i^t 
(27r)2 J J \^{s-i6y + KH^ + d 



X'k 



■dtds. 



The inner integral is [13] 

exp (-i^t 



:dt 



vrexp 



2 



^x^s - i5y + d) 



X'^is - iSy + kH"^ + d K JA2(s - %bf + a 
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and the outer integral is [T3] 

exp f — 7-y/A2(s — i6y + 



X + y 



exp — « — z — s I as 



a'^'Pv 2 



*|A-„|fA'^A-^^ + » 



a; + y 



- exp ^^-5 j 7^0^^ 



/ a 



K^ {x + y) + A^ {x — y) j 



Multiplying by {\ k /{2ti) )) ■ (tt/k), we obtain that f2{x,y) is 



Ak fx + y . 
—- exp ( ^T^d ) Kq 



f ^V «^ (a^ + y? + A2 (x - y)^ 



27r " V 2 

For example, if p = cr = 1/2, the probability that both ^a,c > and 75^(7 > is 
0.608173447.... 



We wish to evaluate 

fsix^y) 

The inner integral is 

exp {-i^t 



72^ 



3. Sample Size n = 3 

exp {-t^s - i^t 



:dt ds. 



(A2(S- 25)2 + ^2^2 + ^^2)3/2 



dt 



(A2(s-i5)2 + fi;2t2 + ^2)3/2 

exp {-i^t 



d 



X^{s-i6)dsJ^ ^A2(S - Z5)2 + fi:2t2 + 



=(it 



a^ 



1 ^2^^ /I 

A2(s — 2(5) (is K Vk 



-X + y 



A2(s-i5)2 + a^ 



K^ 



-X + y 



^^(^ 



2 



A2(s-z5)2 + a2 



^A2(s-z5)2 + a2 



and the outer integral is [13 



^i{7\/A2(s-Z(5)2 + a2j / x + y , ^ 

exp — « — z — s I as 



vr 



^A2(s - i(5)2 + a2 
'x + y 



a;7A 

TT 

a7A 



exp 



exp 



2 

X + y 



a' 



(5 ) exp I —-rxH'^X'^ + 



X + y 



5 ) exp 



a 



2\k 



K^ {x + y) + A2 (x — y)' 



Two Sample Covariances from a Trivariate Normal Distribution 7 



Multiplying by (V2A^/t^/(27r)^) • (27/k), we obtain that fz{x,y) is 
A^K^ fx + y 



1 



exp 



b\ exp f -^V '^^ (3^ + vf + A^ (x - yf\ 



2-1x^-0 



_ exp I - ( px + py - v^V (1 - p2)x2 - 2(-p2 + cj)xy + (1 - p^)!/^ 



which is remarkably simple. For example, if p = a = 1/2, the probability that both 
7a,c > and ^3,0 > is 0.6837762984.... 



We wish to evaluate 



h{x,y) 



The inner integral is 



4. Sample Size n = 4 



2A^4 - - exp(-^2^s-z^t 
(27r)2 



(A2(S - Z5)2 + /s:2t2 + ct2)^ 



rdtds. 



I 



exp f- 



'X+y 



{X^{s-i6y + KH^ + a^] 



idt 




2 



(A2(s-z5)2 + a2)3/2 A2(s-i5)2 + 



and the outer integral is [13 

1 



+ 



.7(A2(s-i5)2 + a2)3/2 A2(s-z5)2 + 



a^ 



exp ( — 7-i/A2(s — Z(5)2 + a"^ \ exp ( —i — - — s ) ds 



a7A2 
1 



exp 



exp 



X + y 
2 

X + y 



2 y^MAV^^^^+V 2 



5j Vfi;2 (x + y)^ + A2 (x - |/)^iri (-jj-\Ji^'^ {x + ?/)^ + A2 (x - y)' 



a'~)X^K '' \ 2 
Multiplying by (2A^K^/(27r)^) ■ (7r7/(2fi;)), we obtain that fA{x,y) is 



Ana 



exp 



X + y 



5) ^«:2(a; + y)2 + A2(x-y)2iri i^^^KHx + yf + XHx~yy 



For example, if p = a = 1/2, the probability that both ■ja,c > and 7b,c > is 
0.7409625593.... 
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5. Formula for Arbitrary n 



We wish to evaluate 



fn{x,y) 



2(n-2)/2^n^n 



(2vr)2 
The inner integral is p3| [T5l [T6] 



exp ^-i^s - i^t 



dtds. 



exp (-«^t) 



dt 



^ 



y (A2(s-25)2 + «:2t2 + «2)n/2 r(n/2)2("-2)/2 «:(n+l)/2 



-X + y 



{n-l)/2 



(A2(s-z5)2 + «2)("-l)/4 



■^(n-l)/2 ( - 



-X + y 



A2(s - iSy + a^ 



and the outer integral is [13 



i^(n-l)/2(7v/A2(g-^^)2 + a2) 

(A2(s-z5)2 + a2)(«-i)/4 



exp —2 



.X + y 



s I (is 



2 a("-2)/2-y(«-l)/2A"/2 

\/7^A2 + ( 



exp 



(i±l,)L.,.,(£±l)^ 



(n-2)/4 



X + y 



(n-2)/2 



exp 



2 2("-4)/2Q,(n-2)/2^(n-l)/2_>^n/2^(n-2)/2 
■^(n-2)/2 ( ^V'^Ma^ + y) +A2(X-?/)' 



X + y 



(j) (/.^ (x + yf + X\x- yf) 



(n-2)/4 



Multiplying by 



2{n~2)/2yn^n 



2n 



(27r)2 
we obtain that fn{x,y) is 

1 A"/2fi:"/2 



:^("-i)/2 



r(n/2) 2"/27ra("-2)A 



exp 



r(n/2)2("-2)/2 fi: 



■^{n-2)/2 f ^V '^^ (^ + 2/)^ + A2 (X - y)^j 



which indeed generalizes the cases n = 1, 2, 3, 4 worked earlier. This result is asymp- 
totically consistent with the following outcome of the Central Limit Theorem: 

i'fA,c-np)/^\^((0\ ( p' + l p' + cr\\ 



ijB,c-np)/y/n I 







' .2 



p^ + a p^ + 1 



as n — > oo. 
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6. Closing Words 
Testing the hypothesis Hq : ■yA,c = 1b,c can be done by examining the difference 



If Hq is true, then 



7A,C - lB,C = 1A-B,C = Y^{.Aj - Bj)Cj. 



'c"h- : ■p-o"? 



and is independent of p; further, the density of ^a-b,c is 

via known results [7] on Gaussian inner products. A considerable literature exists on 
the harder problem of testing Hq : pA,c = Pb,c where variances are unknown and 
underlying distributions might not be normal [HlIISlIISlEniEIlESlESlISlESlESlEZl 

ESlESlEniEIlEaESlEllEaiMlEZlEaiMllinilillliaiia 

We have not attempted to evaluate the inverse Fourier transform of 



G{u, v) = det 




u/2 v/2 ^ ^ ^ '^^^ 
u/2 w/2 
v/2 w/2 



There is no analog of Miller's result [T], as far as we know, giving a joint density 
g{x,y,z) for {AB,AC,BC). Hence no distributional insight on {^A,B,lA,CylB,c) 
seems to be available here. Interestingly, a formula for a joint density for {pa,b, Pa,c, Pb,c) 
is outlined in [511 [52] - evidently a sample size > 4 is presumed - and details still 
need to come together. 
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